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The space-time conservation element solution element (CESE) method is modified to address the robustness 
issues of high-aspect-ratio, viscous, near-wall meshes. In this new approach, the dependent variable gradients 
are evaluated using element edges and the corresponding neighboring solution elements while keeping the 
original flux integration procedure intact. As such, the excellent flux conservation property is retained and 
the new edge-based gradients evaluation significantly improves the robustness for high-aspect ratio meshes 
frequently encountered in three-dimensional, Navier-Stokes calculations. The order of accuracy of the 
proposed method is demonstrated for oblique acoustic wave propagation, shock-wave interaction, and 
hypersonic flows over a blunt body. The confirmed second-order convergence along with the enhanced 
robustness in handling hypersonic blunt body How calculations makes the proposed approach a very 
competitive CFD framework for 3D Navier-Stokes simulations. 

I. Introduction 

I n recent years, the space-time conservation element solution element (CESE) method proposed by Chang [1] has 
attracted cross-discipline attention in computational aerosciences. Applications in different disciplines such as 
complex shock structure flow [2], aeroacoustic wave simulations [3], chemically reacting flow [4], 
magnetohodrodynamics [5], stress waves in solids[6], and the melting process [7-8] have emerged in the past 
decade. The method offers a genuine, multi-dimensional numerical framework for conservation laws that is 
second-order accurate in both space and time. Discretized equations are derived from reinforcing flux conservation 
locally as well as globally using the strong form of the governing equations. Numerical integration of the flux is 
carried out through a set of conservation elements (CE) that do not coincide with the solution elements (SE) where 
the dependent variables are assumed to vary linearly with respect to time and all spatial directions. This staggered 
arrangement of CE and SE completely eliminates the need to solve a Riemann problem along the element interface - 
a practice that has been the central notion of many finite -volume [9-11], spectral difference [12], and discontinuous 
Galerkin (DG) methods [13-15]. The CESE method has been shown to provide remarkable numerical accuracy in 
the cross-discipline application mentioned above. 

Despite decades of successfully applying modern upwind schemes in solving many supersonic problems with 
discontinuities using structured grids, it is becoming apparent that a straight forward extension of the numerical 
framework to calculations using an unstructured mesh may not work in many cases involving strong shocks [16-17], 
In some applications, careful tuning of the flux limiters and entropy fixes may enhance the robustness of the 
schemes for unstructured meshes [18-19]. However, no universal fix is currently available. The problem rests on 
two main issues. The first is the lack of a multi-dimensional extension of the Riemann solver that is generally 
applicable to arbitrary triangular or tetrahedral meshes for the reconstruction of the interface flux. While 
dimensional splitting is often used to make the reconstruction locally one-dimensional for structured meshes, no 
proper counterpart is available for multi-dimensional unstructured meshes. The second issue is associated with the 
famous carbuncle phenomenon appearing in many hypersonic blunt body calculations [20-21]. Numerical 
instability near the bow shock of a blunt body causes the shock around the stagnation region to change shape, and 
the solution eventually diverges. Special treatment in the upwind formulation is required to maintain stability. In 
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some cases, the stabilized numerical solution still exhibits abnormal behavior inside the boundary layer near the 
stagnation region [16]. 

This paper focuses on applying the space-time CESE method for 3D Navier-Stokes calculations. In the past, the 
CESE method has been shown to provide accurate solutions of Navier-Stokes equations for several steady-state or 
time-accurate 2D problems [22], Despite the remarkable success in solving the Euler equations, it has also been 
shown that when high aspect ratio viscous meshes are encountered, the CESE method becomes unstable owing to 
numerical problems with the dependent variable gradients in the elongated element [22], The stiffness with high 
aspect ratio meshes is common to many numerical methods. However, the CESE method appears to be more 
sensitive because the dependent variable gradient plays an important role in the linearly varying solution element. 
The capability to handle a mesh with aspect ratio as high as 10 5 -10 6 is crucial for 3D, Reynolds Averaged Navier- 
Stokes (RANS) or large eddy simulation (LES) calculations because, in general, a mesh with a y + value of 1 or 
smaller is frequently encountered in these simulations. 

The main objective of this paper is to investigate the validity of a modified CESE method proposed to deal with the 
high aspect ratio mesh problems in 3D, Navier-Stokes calculations. The modified scheme is aimed at improving the 
dependent variable gradient calculations that are problematic in high aspect ratio meshes. We discuss the modified 
scheme in detail in the next section. The order of accuracy of the proposed scheme is demonstrated for acoustic 
wave propagation, shock and wave interaction, and bow shocks over blunt bodies in the results section. Finally, 
some conclusions are drawn at the end of the paper. 


II. Numerical Method 

The Original Space-time CESE method 

The numerical formulation of the space-time CESE method has been discussed in detail in the literature [1, 23,-24]. 
We briefly review the numerical formulation below before discussing the modified scheme. The three-dimensional 
compressible Navier-Stokes equations in vector form can be written as 

dQ c)F })G dH „ c)F dG dH v 

— + — + — + = K + — - + — 11 + - (1) 

dt dx dy dz dx dy dz 

where the dependent variable vector is defined as Q = (p, pu, pv, pw, e) and x, y, z, and t represent spatial 
coordinates and time, respectively. Flow variables p, u, v, w, and e represent density, three velocity components, and 

total energy ( e = p l(y — 1) + p(u 2 + V” + W 2 ) / 2 ), respectively. Definitions of the inviscid flux vectors F, G, 
H and viscous flux vectors F„ G v , H v can be found in standard text books (e.g. Hirsch [25]) and will not be included 
here. The source vector K contains all external forcing or other energy-related source terms. For axisymmetric flows, 
H = H v = 0 and K contains extra terms associated with the cylindrical coordinate system. The Euler equations are 
recovered by neglecting all viscous and source terms on the right hand side of the equation. To close the system, the 

perfect gas relation, p — pRT is used in conjunction with eq. (1). 


The space-time CESE method is formulated in the strong form of the flux equations. The space-time flux h is 
defined as: 


h = (F — F V ,G — G V ,H — H V ,Q) (2) 

By using Gauss’ divergence theorem in the space-time domain, eq. (1) is rewritten in the following integral form 

£ h ■ ds + f K dV = 0 (3) 

Js(V) J 
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where the space-time flux vector is integrated over the surface S of an arbitrary space-time volume V. The surface 
normal vector is defined by S = n dA , where dA is the area increment on S, and h is the outward unit normal 
vector. Equation (3) is quite general, and in fact, any conservation law with or without a source vector can be cast 
into this form. Thus, the numerical algorithm devised for eq. (3) can be easily extended for other physical problems. 

For an arbitrary 3D element with a solution point at ( xq , yo, Zo, to), the dependent variable vector is assumed to vary 
linearly within the element: 

Q(x, y, z,t)=Q 0 +Q r (t-t 0 ) + Q x (x -x 0 ) + Q y ( y - y 0 ) + Q z (z - z 0 ) (4) 

where Q 0 is the solution vector at the solution point. The vectors Q t , Q x , Q y , and Q z are analogous to the derivatives 
of the dependent variables. For a higher-order extension, quadratic or third order terms can be added. Theoretically, 
the solution point can be any point within the element. In conventional finite-volume based methods, the solution 
point is usually at the geometric center of the element; however, in the CESE method, the solution point is taken to 
be the geometric center of all surrounding integration volumes called conservation elements. As will be shown later, 
this definition leads to much simpler flux integration equations. The coefficients of the polynomial, Q, t Q x , Q y , and 
Q z , are part of the solutions and need to be determined within each element. Chang points out that the derivatives of 
dependent variables must be treated as unknowns in order to preserve the space-time inversion invariant property 
[26]. The use of eq. (4) results in a formally second-order accurate scheme. 




Figure 1. Conservation elements for a 2D triangular mesh: (a) top view (b) space 
time CE volume 


Without loss of generality, we discuss the numerical formulation using a 2D triangular mesh depicted in Figure 1(a). 
Extension to tetrahedral and other types of elements is straightforward. Point O with coordinates (x 0 , y () , t {] ) is the 

solution point of the element formed by A FBD. For a cell-centered, finite-volume or discontinuous Galerkin 
method, numerical integration is done for A FBD directly. As a result, numerical fluxes must be evaluated at all the 

element interfaces FB, BD, and DF . Unfortunately, numerical fluxes cannot be uniquely determined at these 
element interfaces due to its association with more than one element. For incompressible flows, a simple average of 
dependent variables from both sides of the interface is used to determine the Q and flux vectors accordingly. For 
flows with discontinuities, an approximate Riemann solver is usually applied at the interface. While such Riemann 
solvers generally work well for ID equations, no genuine multi-dimensional Riemann solver is available in the 
literature. A trivial extension using a ID Riemann solver is usually imposed at the element interfaces. In contrast, 
the CESE method uses a different set of integration volumes called conservation elements (CE). Let points G, A, C, 
and E denote the geometric centers of element A FBD and its three surrounding elements, respectively. Numerical 
results, as will be discussed later, indicate that using the in-centers (the intersection of three bi-sectors of three 
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internal angles) [34] in lieu of geometric centers leads to better numerical accuracy. For the triangle A FBD, three 
conservation elements formed by quadrilaterals GBAF, GFED, and GDCB are used for numerical integration of eq. 
(3). In the space-time domain, these quadrilaterals form three rectangular prisms with a height of At / 2 as shown 
in Fig. 1(b). Note that the bottom and top faces of these prisms are located at the old and new time level, 
respectively. 

Two main steps are necessary to evaluate the unknowns of eq. (4) at the new time level. The first step involves 
evaluation of Q r , at the new time level. Equation (4) is valid in the vicinity of the solution point O. The region of 

validity spans over three line segments: GF , GD , and GB at the new time level, and it extends At / 2 above and 
below point O. The flux integral (the first term in eq. (3)) is then carried out for all three CE prisms shown in Fig. 
1(b). By summing up these three flux integrals, we obtain a discretized equation for the triangular element of 
interest. This discretized equation involves known values at the old-time level obtained from integrations over the 
outer-side ( FAA F', ABB A' , etc) and bottom faces ( A'B'G'F ' , G'B'C'D' , etc) of all three prisms. The flux 
integrals of the inner-side faces (GFF'G' ,GBB G' , etc) cancel each other, and the only contribution at the new 
time level is associated with the top faces (GBAF, GFED , and GDCB). If the solution point O is defined as the 
geometric center of the polygon ABCDEF , the discretized flux equation for all three prisms can be simplified to 

(A + A 2 + A 3 )Q 0 = /, + / 2 + I 3 + K v (5) 

where A h A 2 , and A 3 are the surface areas of GBAF, GFED, and GDCB, respectively. Discretized integrals I h I 2 , 
and / 3 are flux integrals for outer-side and bottom faces for three surrounding elements. Source terms in eq. (3) can 
be integrated to form K v . Semi-implicit treatment of this source term is usually done to improve solution accuracy 
and avoid stiffness issues in reacting flows [4], Equation (5) is an algebraic equation and can be solved easily 
without any matrix conversion. 

The second step involves evaluation of Q t , Q x , Q y , and Q z In the non-dissipative a-scheme, these derivatives are 
computed by solving the individual flux equation associated with each CE prism in conjunction with the following 
equation: 


Q t +AQ x +BQ y =K 


where A and B are the Jacobian matrices of the flux vectors F and G, respectively. The above equation can be easily 
derived from eq. (1). Thus, there are a total of four equations for four unknowns. Unfortunately, the a-scheme is 
only neutrally stable for the nonlinear Euler and Navier-Stokes equations and cannot be used in general. 
Nevertheless, the existence of such a dissipation-free scheme serves as a baseline to gauge how much numerical 
dissipation is being added in the c-scheme discussed next. To add dissipation for any nonlinear conservation law, 
derivatives are evaluated using finite-differences. For the triangular element of interest, let 
O x (,Vj , y ! ), 0 2 (x n , y 0 ) and 0 3 ( A', , y 3 ) denote solution points of all three surrounding elements (see Fig. 2). On a 
two-dimensional plane, derivatives in x and y can be uniquely determined if any dependent variable is known at 
three distinct points that are not co-linear. From step one, dependent variable Q 0 is already known. Let Q‘ 0 ( i = 1, 

3) denote dependent variables at neighboring solution points evaluated by using the corresponding solution elements 
at the old time level 


Q‘=Q’+Q;At/2 


(6) 


where Q' (j and Q‘ are known solutions at the old time level. Three sets of derivatives can be obtained 


(Ql,Ql)=fAQ 0 ,Qo,Q 2 o,o, o x ,o 2 ) 

(QlQ 2 y )=fAQo, Ql Qo » o, o 2 ,o 3 ) ( 7 ) 
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(QlQl) = fAQ 0 ,Qo, 05 , 0 , 03 . 4 ) 


The function f d denotes derivative evaluation using the three dependent variables and three solution points given. 

In the c-scheme, the derivatives at the new time level is then determined by applying a weighted average over the 
above three sets of derivatives: 


Q x =W(Ql, Ql Ql ) 

( 8 ) 

Q y =W{Q\, Q;. (Ql) 

( 9 ) 


Details of the weighted averaging procedure can be found in ref.[l, 23]. This weighting function returns the 
equivalent of a simple average for smooth flows and a biased value toward smaller derivatives when large gradients 
are encountered. It is analogous to the flux limiters used in conventional upwind schemes to contain large gradient 
or flow discontinuities. 



Figure 2. Triangular elements along with Figure 3. Triangular elements showing 
neighboring solution points used for derivative neighboring solution points along with zero- 
calculations. dissipation locations used for the CNI scheme. 


Despite the use of a similar weighting function, the CESE method outlined above is fundamentally different from 
the conventional upwind schemes due to the following properties: 

1) Discretizations in space and time are treated equally. This guarantees a uniform order of temporal as well as 
spatial accuracy. 

2) Flux conservation is enforced for all three surrounding CE’s of the triangular element FBD in the above analysis. 
When the discretized flux conservation equation (eq. (5)) is applied repeatedly for all the elements in the domain, 
flux conservation is guaranteed for the entire space-time domain of interest. 

3) Numerical integration of the conservation laws is carried out through staggered conservation elements. The 
interfaces of these conservation elements do not cut through different solution elements where the first-order 
approximation using eq. (4) is valid. As a result, there is no need to formulate a Riemann problem and reconstruct 
fluxes at the interfaces, as is done for almost all finite-volume or DG methods. Aside from the difficulties in 
constructing a “genuine” multi-dimensional Riemann solution, any special interface treatment during flux 
integration produces inconsistency in numerical integration over the entire domain. 
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4) Combining properties 2 and 3, the CESE method offers, arguably, the only genuine multi-dimensional scheme to- 
date based on finite-volume integration using the strong form. As will be demonstrated later, a genuine multi- 
dimensional scheme ensures consistency and accuracy in dealing with strong flow discontinuities. 

5) The CESE method offers a discretized scheme that preserves flux conservation not only in all spatial directions, 
but also in time. It has been shown that the space-time conservation property leads to a much more forgiving 
boundary treatment (e.g. ref.[22]). In many unsteady simulations, the boundary reflection is relatively minor and 
confined within a small region near the boundary. A direct benefit is that the proximity of the truncated boundary is 
not an issue for many practical applications. 


The Courant Number Insensitive (CNI) scheme 

Despite the effectiveness of the c-scheme, it can be shown that numerical dissipation increases as the local CFL 
number decreases [1], For highly non-uniform meshes, numerical accuracy is significantly impacted by large 
numerical dissipation. Chang [27] introduced the Courant number insensitive scheme (CNI) scheme to deal with 
small CFL numbers. In the CNI scheme, the neighboring solution points (O/, O?, and Oi ) used for finite difference 
derivative evaluation are parametrically shifted to control numerical dissipation. Referring to Fig. 3, let points 
Pj , P 1 , and P 3 denote the geometric centers of quadrilaterals GBAF, GFED, and GDC, respectively. Points N t , as 
defined by 

N i =P i +a(O i -P i ) ( 10 ) 

are used for derivative evaluation instead of solution points (9,. The parameter o is used to control numerical 
dissipation by moving the point Nj along the ray P j O i . It can be shown that the points P, correspond to “zero” 

dissipation when the time step is approaching zero [27]. Therefore, by using a a value smaller than unity, numerical 
dissipation can be reduced. On the other hand, a a value of unity is equivalent to the original c-scheme. In his 
investigation of Navier-Stokes solutions, Chang [22] pointed out that using a o value greater than unity will increase 
the numerical dissipation beyond what the c-scheme offers. It is crucial to use a larger a value in a high-aspect ratio 
mesh near viscous wall to maintain numerical instability. Numerical experiments performed by Chang [22] also 
indicate that the use of geometric centers for P, may not be as robust as a scheme that uses the midpoints between 
solution points, i.e., 


P i =(0 + 0 i )/ 2 ( 11 ) 

All solutions presented in ref. [22] were carried out by using Eq. (11), not the geometric centers of surrounding 
conservation elements. 

In the CNI scheme, evaluation of neighboring solutions are now replaced by 

Qo=Qo+ Qi A t / 2 +Q[ (x Nt - x t ) +Q\ (y Ni - y , ) (12) 

where Q x and Q' y are solution derivatives at the neighboring solution element and (x N , y N ) denotes the 
coordinates of the point (V, . Neighboring derivatives can then be evaluated by 

(QlQ\)=f d (Q„QlQlo, n v n 2 ) 

0 Q 2 x ,Q 2 y )=fAQ 0 , Ql Qo’O, n 2 ,n 3 ) (13) 

(QlQl)=fAQo> Ql Qlo, n 3 ,n,) 


The same weighting functions, Eqs. (8-9), can be used to compute the solution derivatives as before. 

For a relatively uniform mesh, the CNI scheme can be used to control numerical dissipation when the CFL number 
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is approaching zero. A o value of 0.01 or even smaller can be used when the CFL number is on the order of 0.01 or 
smaller. However, for a highly stretched mesh with a very high aspect ratio near viscous wall, a small a value in the 
large size mesh region may result in numerical instability. In practice, the a value may not be smaller than 0.3-0. 4. 
An optimal distribution of a with respect to the local CFL number has yet to be determined. Depending upon the 
aspect ratio variation, typical values of o range from on the order of 100 to a minimum value of about 0.4 for 
Navier-Stokes calculations [22]. 

An alternative to the CNI scheme for steady-state calculations is to use local time-stepping. In these calculations, 
the local time step varies from element to element, but the local CFL number is kept at a constant value. While this 
causes some inconsistency in the temporal derivatives, the numerical solution should be the same once a steady-state 
is reached. Local time-stepping can speed up convergence to steady state in most applications. For a highly non- 
uniform mesh, a much larger value of a is needed only for very high aspect ratio elements and should be kept small 
for other elements even if local time-stepping is used. In this case, numerical experiments indicate that by varying a 
according to the mesh size and the aspect ratio (using larger values for larger aspect ratios) can also improve the 
solution accuracy and convergence to steady-state. 

The Edge-based Derivative (EBP) Approach 

The advent of the CNI scheme makes CESE computations on a stretched mesh more accurate if no local time- 
stepping is used. However, numerical problems still exist when a high aspect ratio mesh is encountered. For a 
typical viscous flow calculation, when the aspect ratio near the no-slip wall is more than 100, a large o value on the 
order of 10-100 is needed to stabilize the time-marching procedure. If the aspect ratio goes above 10 4 or so, the time 
iteration becomes unstable no matter how large the o is. This numerical instability is associated with the derivative 
calculation in the second step (eq. (7)), not the flux integration in eq. (5). The situation is better illustrated in Fig. 4 
where a typical viscous mesh is shown. As can be observed, some of the neighboring solution points form an ill- 
shape triangle (e.g. AOOj0 2 and A OOjOj) in conjunction with the element solution point O. Two major problems 

arise due to these ill-shape triangles. First, dependent variable derivatives (e.g. f d (Q 0 , Qq, Qq, 0, O l ,0 2 ) in 

eq. (7)) using these solution points may result in significant error if the three points involved are almost co-linear. 
This causes the derivatives computed in the second stage to have a very large values. Large derivatives in turn cause 
large error in the flux integration for the next time step, and the iteration eventually becomes unstable. 




Figure 4. Conservation elements and 
neighboring solution points for a high aspect 
ratio mesh, showing skewness of the triangles 
used for derivative calculations. 


Figure 5. Conservation elements and 
neighboring points used for derivative 
calculations in the edge-based derivative 
approach. 


7 

American Institute of Aeronautics and Astronautics 


The second, and probably more important issue is related to the shape and location of the triangles used for the 
derivative calculations. As shown in Fig. 4, dependent variable derivatives computed at the current time level would 
be used to calculate fluxes on the interfaces of three quadrilateral-shaped CE’s (depicted with dash lines in the 
figure) at the next time level. For instance, flux integration over the interface on the line segment GB would need to 
use eq. (4) and the computed derivatives, Q x and Q y . However, these derivatives are computed using three triangles 
formed by O 1 O 2 0 3 . None of these three triangles cover the line segment GB entirely. A similar observation can be 
made for line GD. Only the line segment GF is located mostly within the triangles used for derivative calculations. 
For a small aspect ratio mesh, as shown in Fig. 2, this inconsistent shape and location may not have major effect on 
the accuracy of derivatives. However, in a typical viscous calculation, the situation is much worse than that shown in 
Fig. 4. Consequently, inaccurate derivatives used for flux integration may lead to numerical problems during time 
marching. 

A better strategy to improve the accuracy is to use a consistent shape and location for both the derivative evaluation 
and flux vector integration. This notion forms the basis for the modified scheme using edge-based derivatives. In 
this new approach (see Fig. 5), we use the vertices of the triangular element instead of the neighboring solution 
points for derivative calculations. In other words, points 0 / 0 2 , and Ot now coincide with B, F, and D, respectively. 

Equivalently, derivative evaluation is now based on three edges of the triangular element. For the side DB (or 

0 ] 0 2 ) that is located within the neighboring solution element I (with solutions defined by Q ,[ , Q [ , Q' ( , and 


Q 1 at a solution point (x : ,y ; ) ), solutions at the two vertices are given by 

Q) = Go + G/ At / 2 +Q[ (x 0i - x ; ) +Q[ (y 0| - y, ) (14) 

G, 2 = Go + G/ At / 2 +Q[ (x 0i - X , ) +Q[ (y 0i - y, ) (15) 

The corresponding derivatives are then 

(G.:,G.:)=/ rf (Go, Q], Qf,0, o x ,o 2 ) ( 16 ) 

Repeating the above procedure for the remaining two sides with neighboring elements II and III gives, 

(Q" ,Qy)=fAQo’ Ql Ql,0, o 2 ,o 3 ) (17) 

(Qf ,Qf) =f d (Qo’ Qim QnnO, 0 3 ,OQ ( 18 ) 


Finally, the derivatives in the element can be computed by applying a weighted average over the above three sets of 
derivatives just obtained. In contrast to the original CESE method, the current edge-based derivatives approach 
calculates dependent variable derivatives based on the element edges. These edge derivatives allow neighboring 
elements to influence the solution at the derivative level. The main advantage of this new approach is that 
derivatives are evaluated using a shape and location that covers all interfaces of the CE. Comparing Figs. 4 and 5, it 
is evident that the misalignment of the triangles for derivative evaluation and flux integration does not occur for the 
new approach. 

In fact, the use of the side faces and element vertices in the derivative evaluation as outlined above corresponds to a 
zero dissipation state when the time step is approaching zero, analogous to the use of geometric centers of the CE’s 
in the CNI scheme (P, in Fig. 3). Numerical experiments also confirmed this observation. For practical Euler or 
Navier-Stokes calculations, numerical dissipation may be added in a fashion similar to the CNI scheme. Instead of 
using vertices of the triangular element, an adjustable parameter c is used to control the points IV, along the rays 

OO x , OO 1 , and 00, according to 

AT. =0 +<7(0, -O ) (19) 
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Unlike the original CNI scheme discussed above, a o value of unity (using points O, for N, ) now corresponds to 
“zero” numerical dissipation. As shown in Fig. 6, the use of a c value greater than 1 effectively expands the triangles 
used for derivative evaluations. An equivalent CNI scheme for the EBD approach can thus be constructed by 
adjusting the o value according to the local CFL number. 


Validation of the EBD Courant Number Insensitive Scheme 


To test the validity of the Courant number insensitive scheme based on the proposed EBD approach, we compute a 
supersonic flow over a 5° ramp with a stretched mesh as shown in Fig. 7. The solution computed by using the 
original c-scheme with a CFL number of 1 is shown in Fig. 8 as the baseline solution. Due to the disparity in grid 
size, the computed shock near the outflow boundary is significantly smeared. Figure 9 shows the solution computed 
by using the EBD Courant number insensitive method. The shock resolution is improved near the outflow 
boundary. The solution can be further improved by using the local time-stepping (constant CFL) method coupled 
with a locally adjusted rr distribution (see Fig. 10). These results indicate that the concept of CNI scheme still works 
well for the current EBD scheme. A a value of in the range of 1 to 10 is normally used in typical calculations using 
the EBD approach. A value close to 1 while maintaining numerical stability should be used to reduce numerical 
dissipation. It is interesting to note that an asymptotic behavior is observed for large o value for some test cases 
when the local time-stepping procedure is used for a very high aspect ratio mesh. Numerical dissipation appears to 
be unchanged beyond certain value of a. 

The use of eq. (19) also plays an important role in Navier-Stokes calculations. For very high aspect ratio meshes, 
the EBD scheme is much more robust than the original scheme due to the reasons mentioned above. For two 
dimensional calculations, the EBD scheme can handle an aspect ratio as large as 10 5 -10 6 . By using a large o value 
on the order of 100 or larger, numerical stability can be maintained for high-aspect-ratio (on the order of 10 5 ) 
tetrahedral meshes. Overall, the use of EBD modification can significantly improve the robustness of three- 
dimensional Navier-Stokes calculations, especially for tetrahedral meshes. Furthermore, the enhanced stability is 
obtained without any loss of accuracy. 



Figure 6. Construction of Courant number Figure 7. Unstructured triangular mesh for a 

insensitive scheme for the edge-based derivative Mach 3 flow over a 5° ramp 

approach. 
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Figure 8. Density contours computed by 
using the original c-scheme with CFL = 1 . 


Figure 9. Density contours computed by using 
the modified CNI scheme with CFL = 1 . 







Figure 10. Density contours computed by the 
modified scheme with local time stepping 
and locally varying o. 


Figure 11. Oblique acoustic wave 
propagates diagonally through a square 
domain, showing density contours. 


III. Results and Discussion 

Several 2D and 3D test problems were investigated to validate the current EBD scheme. We focus on two main 
issues concerning the application of this modified CESE method: order of accuracy and hypersonic flow 
applications. The order of solution accuracy for the CESE method based on numerical results has yet to be 
determined in the literature. Here, we benchmark the solution accuracy for flows with and without shocks; and for 
steady as well as time-accurate calculations. All CESE solutions presented below were computed by using the 
modified EBD scheme. 
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A. Order of Accuracy Study 


To benchmark the order of accuracy for time-accurate calculations, we first compute an oblique acoustic wave 
propagating through a square domain as shown in Fig. 11. No mean flow is present in the domain. Acoustic wave 
inflow is specified at the left and lower boundaries and non-reflecting boundary conditions are applied at the upper 
and right boundaries. An unstructured mesh is generated by triangulating a corresponding uniformed structured 
quadrilateral mesh. For all calculations, the CFL number is fixed at 1. Thus, the time step decreases as the grid is 
refined. The density contours after the wave has completely traversed the domain are shown in the figure. The 
computed solutions within a square box centered at the origin with a length of 1 are compared with the exact 
solution in order to calculate the L2 error norm shown in Fig. 12. The quantity N in the figure is the number of 
points in both the x and y directions for the corresponding structured mesh. When the grid size is doubled, the error 
norm decreases by roughly a factor 4 or slightly better. This error norm confirms that the present scheme is indeed 
second order in both space and time. Although not shown here, the original c-scheme produces similar error norm 
convergence to the present EBD scheme. 

A more challenging benchmark problem to test the solution accuracy includes a shock in the flowfield. A Mach 3 
flow over a 5° ramp is computed with a triangular mesh generated by the corresponding structured meshes 41x21, 
81x41, 161x81, and 321x161. The resulting area-averaged LI error norm is plotted against the number of grid 
points in the stream wise direction in Fig. 13. The general observation is that the error is the largest near the oblique 
shock and small in the smooth region. However, the area-averaged norm is dominated by a large area of smooth 
flow region. The overall LI -norm for large N is approaching second order. This example indicates that the 
weighted average applied in the vicinity of the shock does not have an adverse effect on the solution accuracy away 
from it. This is a necessary property of any shock capturing scheme. 




Figure 12. Convergence of error norm versus Figure 13. Convergence of error norm 

number of points in the mesh for acoustic wave versus number of grid points for a Mach 3 

propagation shown in Fig. 1 1 . flow over a 5 0 wedge. 


A more rigorous test is the accuracy of a time-accurate calculation in the presence of a shock. We investigate the 
one-dimensional acoustic wave and shock interaction problem normally used to assess the accuracy of the 
essentially non-oscillating (ENO) schemes [28]. In this case, an acoustic wave propagates through a Mach 3 moving 
shock wave and generates waves and shocklets behind the shock. The computed density distribution at the non- 
dimensional time of 1.8 using 1600 points is shown in Fig. 14 along with the ENO-RF-3 solution by Shu[28] using 
the same grid resolution. The agreement between both sets of solutions is good. Both the shorter scale near the 
shock and N-shaped shocklets further upstream are properly resolved by the current second order scheme. Figure 15 
shows the computed LI error norm (computed with the solution obtained by using 3200 points as the reference) 
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versus number of points used. From the rate of change of the curve, the solution is only slightly better than first 
order. Apparently, the shock capturing weighted average results in reduced solution accuracy near the discontinuity. 
This result appears to be in line with Casper and Carpenter’s [29] conclusion drawn from time-accurate wave/shock 
interaction inside a ID nozzle. First-order accuracy in the presence of flow discontinuity may appear to be 
unsatisfactory. However, the fact that only a slightly higher than first-order convergence rate was obtained by the 
fifth-order weighted ENO (WENO) scheme (Pirozzoli [30]) implies that the present result is as good as could be 
expected for this difficult problem. By visually comparing the current solutions with those given in ref. [28], it is 
evident that the present solution using 800 points is comparable with the one obtained with 400 points using the 
ENO-RF-3 scheme, and is much better than the 800 point solution computed with the MUSCL type TVD scheme 
shown in the paper. 



x 

Figure 14. Density distribution at t = 1.8 for the 
one-dimensional acoustic wave and Mach 3 
moving shock interaction. 



N 


Figure 15. LI error norm versus grid resolution 
for the ID acoustic wave/shock interaction 
shown in Fig. 14. 


Another test problem usually used to benchmark the accuracy of the ENO scheme is the 2D shock-vorticity wave 
interaction investigated in refs. [28, 30]. In this test case, a vorticity wave interacts with a Mach 8 moving shock 
wave in a rectangular domain. A shock is initially located at x = - 1 and is moving into a region where the flow state 
is defined by 

P= 1 

u = - ^[y sin 9 cos {kx cos 9 + ky sin 0) 
v = yfy sin 0 cos (kx cos 9 + ky sin 9) 

P= 1 

The values of k and 0 are 2;t and n/6, respectively. Periodic boundary conditions are applied at the top and bottom 
boundaries. A triangulated uniform 384x256 structured mesh is used for computations. Figure 16 shows the 
computed density contours at a non-dimensional time of 0.2. The comparison with the exact solution along the y = 0 
line is depicted in Fig. 17. The wave structure is properly resolved by using the current scheme, and the agreement 
with the exact solution is good except for missing the sharp gradient near x = 0.55. Comparing with solutions 
obtained by the WENO scheme in ref. [30] indicates that with half of their mesh size, the present solution matches 
roughly their W5 (fifth-order) WENO scheme. The 7-th order schemes of [30] give a slightly better solution near x = 
0.55. 
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Figure 16. Computed density contours for the 
shock-vorticity wave interaction problem at t = 
0.2, the shock is initially located at x = -1. 


Figure 17. Comparison of density distribution 
along y = 0 with exact solution for the shock- 
vorticity wave interaction problem shown in 
Fig. 16. 


B. Hypersonic Blunt Body Computations 

In the previous section, we have demonstrated the solution accuracy of the present scheme for several time-accurate 
and steady-state benchmark problems with flow discontinuities. This section is focused on hypersonic, blunt-body 
computations. Accurately computing hypersonic blunt body flows and maintaining numerical stability is one of the 
outstanding issues related to hypersonic flow simulations using unstructured meshes. Despite decades of success in 
computing bow shocks using upwind schemes with structured meshes, there are still reports of the famous carbuncle 
phenomenon in which numerical instability associated with the stagnation region behind the shock causes the shock 
to become unstable [16-17, 20-21], For the configurations where stable solutions can be obtained, the carbuncle 
phenomenon may still be observed if more grid points are used between the shock and body in the calculation [31]. 
For the finite-volume or DG based unstructured methods, this numerical instability appears to be a more serious 
issue (e.g. refs. [16-17]). Due to the scope limitation, we only focus on two issues related to blunt body 
computations: 1) the order property of the present scheme in calculating bow shocks 2) assessment of the solution 
quality for hypersonic flow over a cylinder with strong surface heat transfer. 

A Mach 6 flow over a circular cylinder is chosen to study the order of accuracy of the present scheme. The baseline 
solution used to benchmark the order of accuracy is obtained by a shock-fitting code (Salas and Atkins [32]). In 
such shock-fitting calculations, a discretized domain between the shock and the body is used in solving the Euler 
equations. The main advantage of this approach is that any issues related to numerical instability and the accuracy 
constraint of the shock capturing scheme do not play any role. If sufficient points are used, a “clean” solution can 
thus be regarded as an exact solution. Structured meshes with dimensions 31x41 (31 along the cylinder surface and 
41 along the wall-normal direction), 61x81, 121x161, and 241x321 are triangulated for the present CESE 
calculations. Figure 18 shows the 31x41 mesh near the symmetry line. In these Euler calculations, the slip wall 
condition is used on the cylinder surface. Free-stream and symmetry conditions are applied at the inflow and center 
line, respectively. Figures 19 and 20 compare the computed Mach numbers and pressure contours using the finest 
grid (241x321) with the shock-fitting solution using a 192x320 mesh. Excellent agreement between these two 
methods is evident. 

The shock-fitting solution provides an excellent baseline for the validation of several variations of the CESE 
method. Numerical results indicate that using the in-centers instead of the geometric centers of the triangular 
element to construct the surrounding CE’s (points G, A, C, and E in Fig. 1) is apparently more accurate. A 
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triangular mesh is also more accurate than a comparable quadrilateral mesh with the same number of elements for 
this Mach 6 bow shock problem. Figure 21 depicts the L2 error norm for various mesh sizes. These error norms 
were calculated at the shock-fitting calculation grid points on the cylinder surface and along the symmetry line. The 
CESE solutions at these discrete points were calculated by using eq. (4). As suggested in ref. [32], using an integral 
quantity for the error norm calculation is to be avoided to ensure appropriate measure of order properties. The rate 
of change of both curves shown in Fig. 21 indicates only about first-order accuracy. An area-averaged error norm as 
computed for the supersonic flow over a wedge (Fig. 13) is not available for this case. The reason for the lower than 
second order accuracy is unknown. However, the non-orthogonal mesh near the stagnation region (see Fig. 18) may 
affect the solution accuracy of the present scheme. 



Figure 18. Triangular mesh used for the Mach 6 
bow shock calculation, showing the coarsest 
31x41 mesh. 



Figure 19. Comparison of the Mach number 
contours computed by the shock-fitting code 
and the present method for the Mach 6 flow 
over a cylinder. 



Figure 20. Comparison of the pressure contours 
computed by the shock-fitting code and the 
present method for the Mach 6 flow over a 
cylinder. 


Figure 21. L2 error norm on the cylinder 
surface and the symmetry line versus number 
of points (along the cylinder surface) used for 
the Mach 6 flow over a cylinder case. 


14 

American Institute of Aeronautics and Astronautics 


As pointed out in [16], aerothermodyamic simulation of a blunt body using unstructured grids still poses numerical 
issues associated with the carbuncle phenomenon and the ineffectiveness of the entropy fixes. To accurately predict 
surface heat transfer, a highly refined mesh near the body surface must be used. The configuration used by the 
LAURA benchmark [16] is a Mach 17 flow over a 2D or 3D cylindrical surface. Figure 22(a) shows the triangulated 
61x65 mesh used for simulation. The highest aspect ratio near the wall is about 1400, and, as a result, a o value on 
the order of 100 is needed to achieve a converged solution. Figure 22(b) shows the computed temperature contours 
around the body after the L2 norm has dropped 5 orders of magnitude. As shown, the mesh used is intended to 
properly resolve both the bow shock and thermal boundary layer. Figure 23 shows a comparison of cylinder surface 
heating and non-dimensional pressure between the present results and that from the LAURA code [16]. The overall 
agreement is marginal, especially with the surface heating being under-predicted. However, the present solution is 
by no means a grid-converged solution. Another run with a higher grid stretching near the wall while keeping the 
same 65 points in the wall-normal direction appears to over-predict the surface heating rate. This indicates that the 
solution shown here is not grid-converged. Nevertheless, the smooth surface heating and pressure distributions in 
Fig. 23 confirm that the present method does not have numerical stability problems in handling hypersonic flows 
over a 2D blunt body flows with strong surface heat transfer. 



(a) (b) 


Figure 22. Triangular mesh and computed 
temperature contours for a Mach 1 7 flow over a 
cylinder: (a) mesh (b) temperature contours. 



Figure 23. Comparison of surface heating and 
non-dimensional surface pressure with 
LAURA results [16]: symbols for LAURA 
and lines for CESE. 


To demonstrate that the scheme does not suffer from stability problems in three dimensions, the above 2D mesh is 
stacked up 10 times along the spanwise direction to form a tetrahedral 3D mesh as shown in Fig. 24(a). Periodic 
boundary conditions are imposed at the spanwise boundary. The computed temperature contours are shown in Fig. 
24(b). Although not visible on the scale of the figure, a weak spanwise gradient does exist in the solution. The 
largest spanwise gradient is taking place between the first two planes near the spanwise boundaries. Using a 
reflecting (slip wall) boundary condition does not seem to help in eliminating this weak gradient. The directional 
bias of a tetrahedral mesh may have contributed to this spanwise gradient, but more studies are necessary to confirm 
this. Figure 25 shows a comparison of the computed 3D surface heat transfer and non-dimensional pressure with 
those obtained from the 2D simulation shown in Fig. 23. Clearly, the weak spanwise gradient does affect the overall 
solution. In light of the differences in the discretized volumes between the 2D and 3D simulations (unless a 
prismatic mesh is used, the 3D tetrahedral mesh does not replicate the 2D mesh), we do expect some sort of 
discrepancies. Unlike the surface pressure, the heat flux computed at various spanwise planes appears to be varying. 
The middle planes away from the boundaries seem to agree well and peak around 40 W/cnr. This value is even 
smaller than that from the 2D results. To obtain a grid converged solution is beyond the scope of this paper. 
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However, the present 3D simulation using a tetrahedral mesh does not exhibit any numerical issues concerning 
numerical instability nor errant behavior in the surface heating near the stagnation region. As mentioned in the 
previous section, the CESE method does not have any intrinsic flux splitting issues near the stagnation region. The 
weighted average procedure for the dependent variable gradients appears to be free from the glitches of conventional 
upwind schemes. As such, a genuine multi-dimensional formulation coupled with a unified weighted averaging for 
dependent variable gradients in the present CESE method ensures that no improper entropy generation near the bow 
shock is allowed to contaminate the surface heating and pressure. 

Although not shown here, numerical experiments with a realistic triangular or tetrahedral mesh do not exhibit any 
numerical difficulties either. A mesh that is misaligned with the bow shock also can be handled very well. As 
discussed above, the issues associated with accurate surface heating prediction need to be sorted out. Nevertheless, 
the present method does offer a more robust and viable approach to compute hypersonic blunt body configurations. 


z 



x (m) 


Figure 24. Tetrahedral mesh and computed Figure 25. Comparison of computed surface 

temperature contours for a Mach 17 flow heating and non-dimensional surface pressure 

over a cylinder: (a) mesh (b) contours. using a 2D triangular or a 3D tetrahedral mesh. 


C. 3D Navier-Stokes Simulations 

With the modified scheme, the CESE method now can handle very high aspect ratio viscous meshes without 
numerical instability. For a very high aspect ratio mesh beyond 10 5 , a large c value on the order 100-1000 must be 
used. This implies large numerical dissipation. Solution accuracy with such a large o value has not yet been 
assessed. However, many applications with an aspect ratio below 10 5 can be computed with properly controlled 
numerical dissipation. The modified scheme has been fully implemented in the ez4d code [22]. Pending validation, 
both Spalart-Allmaras and k-omega turbulence models are also available in this code. 

Numerical tests have been performed with the modified CESE method for several 3D Navier-Stokes calculations 
using tetrahedral and hexahedral meshes. Figures 26(a) and 26(b) show the tetrahedral mesh and computed pressure 
contours for a test case of transonic flow (with a free-stream Mach number of 0.85 and -6° angle of attack) over a 
wing body configuration that has been used as a benchmark problem for the CFD Drag Prediction Workshop [33]. 
The results were obtained with laminar calculations. Turbulent calculations along with comparison with other drag 
prediction workshop solutions will be performed in the future. 
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Figure 26. Tetrahedral mesh and computed pressure contours on a three-dimensional wing 
body configuration using the modified CESE method and the ez4cl code. 


IV. Conclusion 

A new procedure to calculate the dependent variable gradients for the CESE method has been proposed. Instead of 
using the neighboring solution points, the modified approach uses the neighboring solution elements and the vertices 
of the element itself to evaluate derivatives influenced by each neighbor. The weighted average of all computed 
neighboring derivatives is then imposed as the resulting dependent variable gradient within each element. The shape 
and location used for the derivative calculations in this modified approach are similar to those used for the flux 
integration. As a result, dependent variable gradients may be more accurately computed for high-aspect-ratio 
meshes, and the time advancement procedure becomes much more robust. It has been shown that a comparable 
Courant number insensitive scheme and numerical dissipation control may be established by expanding the element 
size. A larger but similar element brings about more numerical dissipation to the solution. Numerical dissipation 
can be controlled by using a single parameter that determines the size of the element used for derivative calculations. 

Solution accuracy of this new approach has been validated by investigating the error norm convergence for several 
benchmark problems: oblique acoustic wave propagation in a square domain, Mach 3 flow over a wedge, ID shock- 
acoustic wave interaction, and 2D shock-vorticity wave interaction problems. The results indicate that for steady 
state and shock-free time-accurate problems the modified scheme is second-order accurate. For a time-accurate 
wave and shock interaction problem, the scheme is only first order accurate. This reduced order behavior for wave- 
shock interaction is comparable with that observed for a higher order ENO or WENO scheme. 

Hypersonic flow over a cylinder is used to assess capabilities of the proposed method in handling bow shock flows 
robustly and accurately. Results for a Mach 6 flow over a cylinder are shown to compare very well with an “exact” 
solution obtained by a shock-fitting scheme. The solution accuracy along the symmetry line and the cylinder 
surface has been shown to be only of first order accuracy despite the steady-state nature of this test problem. 
Computations of a Mach 17 flow over a cylinder have been performed to assess the capability in surface heating 
prediction. For a 61x65 mesh, the present results are close to those obtained by the structured LAURA code[16]. 
Both 2D triangular and 3D tetrahedral results show consistent shock profiles and smooth surface heating 
distributions. Neither the carbuncle phenomenon nor abnormal surface heating distributions due to numerical issues 
regarding shock capturing and entropy fixes are observed. Implementation of this modified scheme in a 3D Navier- 
Stokes code shows that the present numerical framework provides a more robust method for general hypersonic 
flow simulations. By and large, the established solution accuracy benchmarks for acoustic wave propagation, shock 
and wave interaction, and blunt body flows demonstrates that the present CESE method offers a robust and accurate 
numerical framework for hypersonic computations using unstructured grids. 
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